###############################################################################
# Replication Code for All-Party Analysis 

# "Strengthening mainstream consensus? The effect of radical right parties on the
# defense policies of left parties"

# Miku Matsunaga & Thomas Winzen 

###############################################################################
sink(file = "log_file1.txt", type = "output")

require(tidyverse) 
require(rio)      
require(magrittr)
require(dplyr)
require(stringr)
require(ggplot2)
require(tidyr)
require(MASS)
require(stargazer)
require(coefplot)
require(sjPlot)
require(margins)
require(plm)
require(estimatr)
require(lmtest)
require(lme4)
require(pscl)
require(jtools)
library(broom)     
library(rdrobust)   
library(rddensity)
library(rddtools)
library(huxtable)   
library(rdd)
require(haven)
require(plyr)
require(cowplot)
library(data.table)
require(gridExtra)
library(estimatr)
library(readr)
require(grid)
require(ggpubr)

###############################################################################

#  Main Analysis 

###############################################################################

#download data
manifesto_pop<-read.csv('data.csv')

#descriptive statistics (Table A1)
stargazer(manifesto_pop)

#creating variable for country FE
manifesto_pop$state.f = factor(manifesto_pop$countryname)
manifesto_pop$state.d = model.matrix(~manifesto_pop$state.f+0)

# variables 
#per104 Military Positive
#per105 Military Negative
#per1011 US Positive
#per1012 US Negative 
#per1021 Russia Positive
#per1022 Russia Negative
#military: defense position (per104-per105)
#military_change: change of defense position (between t and t-1)
#military2: defense salience (per104+per105)
#us:salience of united states  
#russia: salience of Russia
#dif: difference between RRPP vote and electoral threshold 
#difl_1 : difference between RRPP vote and electoral threshold(t-1)
#dif_fixed1: difference between RRPP vote and legal electoral threshold(t-1)


#-----------------------------------------------------------------------------
# 1.Position: Military Positive - Military Negative Per104-Per105 (Table 1, Figure 1, Table A4)
#-----------------------------------------------------------------------------
#%%%%%%%%%  Table 1. %%%%%%%%%%%%%%%%%%%
manifesto_pop$military_l<-log(manifesto_pop$military+1)
manifesto_pop1<-manifesto_pop%>%dplyr::select(dif_l1, military_l,state.d, countryname,treatment)
manifesto_pop1<-manifesto_pop1[complete.cases(manifesto_pop1), ]
mil_f<- rdd_data(x=dif_l1, y=military_l, z=manifesto_pop1$treatment, covar=cbind(state.d), cutpoint=0, data=manifesto_pop1) 
mil_f2<- rdd_data(x=dif_l1, y=military_l, z=manifesto_pop1$treatment, cutpoint=0, data=manifesto_pop1) 

#parametric (fuzzy with fixed effect)
summary(f<-rdd_reg_lm(rdd_object=mil_f, covar.opt = list(strategy = c("include"), slope = c( "separate")), order=1))  #fuzzy w/FE
coeftest(f, vcov=vcovHC(f,type="HC0",cluster="countryname"))

#parametric (fuzzy without fixed effect)
summary(f2<-rdd_reg_lm(rdd_object=mil_f2, order=1)) 
coeftest(f2, vcov=vcovHC(f2,type="HC0",cluster="countryname"))

#nonparametric without fixed effect 
summary(rdbwselect_2014(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_l1, c=0, bwselect="CCT"))
summary(rdrobust(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, p=2, all=TRUE, covs=cbind(manifesto_pop1$state.d),h=2.051, b=3.763,cluster=manifesto_pop1$countryname, fuzzy = manifesto_pop1$treatment))
#nonparametric with fixed effect
summary(rdrobust(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, p=2, all=TRUE,h= 2.051, b=3.763, cluster=manifesto_pop1$countryname, fuzzy = manifesto_pop1$treatment))


########% Figure 1. %#########
d1<-manifesto_pop1%>%dplyr::filter(dif_l1>-5, dif_l1<5)
pdf("Fig.1.pdf")
rdplot(x=d1$dif_l1, y=d1$military_l, c=0, p=2, covs=d1$state.d,x.label="Difference between RRPP vote and threshold (t-1)", y.label="Position of National Defence" , title="")
dev.off()
########% Table A4. %#########
#testing with different bandwidth 
#MSE (one common MSE-optimal bandwidth)
summary(rdrobust(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, p=2, all=TRUE, covs=cbind(manifesto_pop1$state.d),bwselect="mserd",cluster=manifesto_pop1$countryname, fuzzy = manifesto_pop1$treatment))
#CER (one common CER-optimal bandwidth selector for the RD treatment effect estimator)
summary(rdrobust(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, p=2, all=TRUE, covs=cbind(manifesto_pop1$state.d),bwselect="cerrd",cluster=manifesto_pop1$countryname, fuzzy = manifesto_pop1$treatment))
#IK
summary(rdrobust(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, p=2, all=TRUE, covs=cbind(manifesto_pop1$state.d),h=4.142, b=4.376,cluster=manifesto_pop1$countryname, fuzzy = manifesto_pop1$treatment))


#-----------------------------------------------------------------------------
# Placebo test --- arbitrarily set 5% cut off point (Table A6)
#-----------------------------------------------------------------------------
manifesto_pop1<-manifesto_pop%>%dplyr::select(dif_l1, military_l,state.d, countryname, treatment_placebo)
manifesto_pop1<-manifesto_pop1[complete.cases(manifesto_pop1), ]
mil_f<- rdd_data(x=dif_l1, y=military_l, z=manifesto_pop1$treatment_placebo, covar=cbind(state.d), cutpoint=5, data=manifesto_pop1) 
mil_f2<- rdd_data(x=dif_l1, y=military_l, z=manifesto_pop1$treatment_placebo, cutpoint=5, data=manifesto_pop1) 

#parametric 
summary(f<-rdd_reg_lm(rdd_object=mil_f, covar.opt = list(strategy = c("include"), slope = c( "separate")), order=1))
coeftest(f, vcov=vcovHC(f,type="HC0",cluster="countryname"))
summary(f2<-rdd_reg_lm(rdd_object=mil_f2, order=1))
coeftest(f2, vcov=vcovHC(f2,type="HC0",cluster="countryname"))#fuzzy w/o FE

#nonparametric
summary(rdbwselect_2014(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_l1, c=5, bwselect="CCT"))
summary(rdrobust(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_l1, c=5, kernel = "tri",level = 95, p=2, all=TRUE,h=4.628, b=6.646,cluster=manifesto_pop1$countryname, fuzzy = manifesto_pop1$treatment_placebo, covs=cbind(manifesto_pop1$state.d)))
#without FE
summary(rdrobust(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_l1, c=5, kernel = "tri",level = 95, p=2, all=TRUE,h=4.628, b=6.646, cluster=manifesto_pop1$countryname, fuzzy = manifesto_pop1$treatment_placebo))

#-----------------------------------------------------------------------------
# Legal (fixed) threshold only (do it with sharp design) (Table A7)
#-----------------------------------------------------------------------------
manifesto_pop1<-manifesto_pop%>%dplyr::select(dif_fixed1, military_l,treatment_fixed, state.d, countryname)
manifesto_pop1<-manifesto_pop1[complete.cases(manifesto_pop1), ]
mil_f<- rdd_data(x=dif_fixed1, y=military_l,covar=cbind(manifesto_pop1$state.d), cutpoint=0, data=manifesto_pop1) 
mil_n<- rdd_data(x=dif_fixed1, y=military_l, cutpoint=0, data=manifesto_pop1) 

#sharp para
summary(f<-rdd_reg_lm(rdd_object=mil_f,covar.opt = list(strategy = c("include"), slope = c( "separate")), order=1 )) #sharp w/FE 
coeftest(f, vcov=vcovHC(f,type="HC0",cluster="countryname"))

summary(n<-rdd_reg_lm(rdd_object=mil_n,order=1 )) #sharp w/o FE
coeftest(n, vcov=vcovHC(n,type="HC0",cluster="countryname"))

#nonparametric sharp
summary(rdbwselect_2014(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_fixed1, c=0, bwselect="IK", all=TRUE)) 
#with FE
summary(rdrobust(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_fixed1, c=0, kernel = "tri",level = 95, p=2,h=1.395,b=2.609, covs=cbind(manifesto_pop1$state.d), cluster=manifesto_pop1$countryname, fuzzy = manifesto_pop1$treatment_fixed, all=TRUE))#fuzzy
#without FE
summary(rdrobust(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_fixed1, c=0, kernel = "tri",level = 95, p=2,h=1.395,b=2.609, cluster=manifesto_pop1$countryname, fuzzy = manifesto_pop1$treatment_fixed, all=TRUE))

#-----------------------------------------------------------------------------
# Different polynomial order (Table A5)
#-----------------------------------------------------------------------------
manifesto_pop1<-manifesto_pop%>%dplyr::select(dif_l1, military_l,state.d, countryname,treatment)
manifesto_pop1<-manifesto_pop1[complete.cases(manifesto_pop1), ]
mil_f<- rdd_data(x=dif_l1, y=military_l, z=manifesto_pop1$treatment, covar=cbind(state.d), cutpoint=0, data=manifesto_pop1) 
mil_f2<- rdd_data(x=dif_l1, y=military_l, z=manifesto_pop1$treatment, cutpoint=0, data=manifesto_pop1) 

#parametric (fuzzy with fixed effect)
# polynomial =1 
summary(d1<-rdd_reg_lm(rdd_object=mil_f, covar.opt = list(strategy = c("include"), slope = c( "separate")), order=1))  #fuzzy w/FE
coeftest(d1, vcov=vcovHC(d1,type="HC0",cluster="countryname"))
#polynomial =2 
summary(d2<-rdd_reg_lm(rdd_object=mil_f, covar.opt = list(strategy = c("include"), slope = c( "separate")), order=2))  #fuzzy w/FE
coeftest(d2, vcov=vcovHC(d2,type="HC0",cluster="countryname"))
#polynomial =3 
summary(d3<-rdd_reg_lm(rdd_object=mil_f, covar.opt = list(strategy = c("include"), slope = c( "separate")), order=3))  #fuzzy w/FE
coeftest(d3, vcov=vcovHC(d3,type="HC0",cluster="countryname"))
#polynomial =4
summary(d4<-rdd_reg_lm(rdd_object=mil_f, covar.opt = list(strategy = c("include"), slope = c( "separate")), order=4))  #fuzzy w/FE
coeftest(d4, vcov=vcovHC(d4,type="HC0",cluster="countryname"))

#nonparametric without fixed effect 
summary(rdbwselect_2014(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_l1, c=0, bwselect="CCT"))
#nonparametric with fixed effect
# polynomial =1
summary(rdrobust(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, p=1, all=TRUE,h= 2.051, b=3.763, cluster=manifesto_pop1$countryname, fuzzy = manifesto_pop1$treatment))
# polynomial =2
summary(rdrobust(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, p=2, all=TRUE,h= 2.051, b=3.763, cluster=manifesto_pop1$countryname, fuzzy = manifesto_pop1$treatment))
# polynomial =3
summary(rdrobust(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, p=3, all=TRUE,h= 2.051, b=3.763, cluster=manifesto_pop1$countryname, fuzzy = manifesto_pop1$treatment))
# polynomial =4
summary(rdrobust(y = manifesto_pop1$military_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, p=4, all=TRUE,h= 2.051, b=3.763, cluster=manifesto_pop1$countryname, fuzzy = manifesto_pop1$treatment))


################################################################################
#-----------------------------------------------------------------------------
# Salience: Military Positive + Military Negative Per104+Per105 (Table A8 )
#-----------------------------------------------------------------------------
manifesto_pop$military2_l<-log(manifesto_pop$military2+1)
manifesto_pop1<-manifesto_pop%>%dplyr::select(dif_l1, military2_l,state.d,countryname, treatment)
manifesto_pop1<-manifesto_pop1[complete.cases(manifesto_pop1), ]
mils_f<- rdd_data(x=dif_l1, y=military2_l,covar=cbind(state.d), z=manifesto_pop1$treatment, cutpoint=0, data=manifesto_pop1) 
mils_f2<- rdd_data(x=dif_l1, y=military2_l, z=manifesto_pop1$treatment, cutpoint=0, data=manifesto_pop1)

#parametric fuzzy
summary(t<-rdd_reg_lm(rdd_object=mils_f,covar.opt = list(strategy = c("include"), slope = c( "separate")), order=1))
coeftest(t, vcov=vcovHC(t,type="HC0",cluster="countryname"))
summary(t2<-rdd_reg_lm(rdd_object=mils_f2, order=1))
coeftest(t2, vcov=vcovHC(t2,type="HC0",cluster="countryname"))


#non-parametric fuzzy
summary(rdbwselect_2014(y = manifesto_pop1$military2_l, x = manifesto_pop1$dif_l1, c=0, bwselect="CCT"))
summary(rdrobust(y = manifesto_pop1$military2_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, p=2,  h=1.811, b=3.085, covs=cbind(manifesto_pop1$state.d), cluster=manifesto_pop1$countryname, fuzzy = manifesto_pop1$treatment, all=TRUE))
summary(rdrobust(y = manifesto_pop1$military2_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, p=2,  h=1.811, b=3.085,  cluster=manifesto_pop1$countryname, fuzzy = manifesto_pop1$treatment, all=TRUE))


################################################################################
#-----------------------------------------------------------------------------
# Salience:  Russian Salience: per1011+per1021”: Russian positive + Russian negative
#-----------------------------------------------------------------------------

#Table A2 

manifesto_pop$russia_l<-log(manifesto_pop$russia+1)
manifesto_pop1<-manifesto_pop%>%dplyr::select(dif_l1, russia_l,state.d, countryname, treatment)
manifesto_pop1<-manifesto_pop1[complete.cases(manifesto_pop1), ]
rus_f<- rdd_data(x=dif_l1, y=as.numeric(russia_l),covar=cbind(state.d), z=manifesto_pop1$treatment, cutpoint=0, data=manifesto_pop1) #fuzzy
rus_n<- rdd_data(x=dif_l1, y=as.numeric(russia_l), z=manifesto_pop1$treatment, cutpoint=0, data=manifesto_pop1) #fuzzy

#parametric fuzzy
summary(r2<-rdd_reg_lm(rdd_object=rus_f,covar.opt = list(strategy = c("include"), slope = c( "separate")), order=1))
coeftest(r2, vcov=vcovHC(r2,type="HC0",cluster="countryname"))
summary(r3<-rdd_reg_lm(rdd_object=rus_n, order=1))
coeftest(r3, vcov=vcovHC(r3,type="HC0",cluster="countryname"))

#nonparametric fuzzy
summary(rdbwselect_2014(y = manifesto_pop1$russia_l, x = manifesto_pop1$dif_l1, c=0, bwselect="CCT"))
summary(rdrobust(y = manifesto_pop1$russia_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, covs=manifesto_pop1$state.d,fuzzy=manifesto_pop1$treatment, p=2, h= 1.328493, b=4.129771, cluster=manifesto_pop1$countryname, all=TRUE))
summary(rdrobust(y = manifesto_pop1$russia_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95,fuzzy=manifesto_pop1$treatment, p=2, h= 9.443, b=5.200, cluster=manifesto_pop1$countryname, all=TRUE))

#plot (Figure A1)
d1<-manifesto_pop1%>%dplyr::filter(dif_l1>-4.477, dif_l1<4.477)
pdf("Fig.A1a.pdf")
rdplot(x=d1$dif_l1, y=as.numeric(d1$russia), covs=d1$state.d, c=0, p=2, x.label="", y.label="Russian Salience" , title="",y.lim = c(-0.1,0.3), x.lim = c(-5,5))
dev.off()

################################################################################
#-----------------------------------------------------------------------------
# US Salience: per1012/1022 measure salience US and other western countries
#manifesto_pop$us<-manifesto_pop$per1012+manifesto_pop$per1022
#-----------------------------------------------------------------------------

#Table A2
manifesto_pop$us_l<-log(manifesto_pop$us+1)
manifesto_pop1<-manifesto_pop%>%dplyr::select(dif_l1, us_l,state.d, countryname, treatment)
manifesto_pop1<-manifesto_pop1[complete.cases(manifesto_pop1), ]
usa_f<- rdd_data(x=dif_l1, y=as.numeric(us_l), z=manifesto_pop1$treatment, covar=cbind(manifesto_pop1$state.d),cutpoint=0, data=manifesto_pop1) 
usa_f2<- rdd_data(x=dif_l1, y=as.numeric(us_l), z=manifesto_pop1$treatment,cutpoint=0, data=manifesto_pop1)


#parametric fuzzy
summary(i<-rdd_reg_lm(rdd_object= usa_f, covar.opt = list(strategy = c("include"), slope = c( "separate")),order=1)) 
coeftest(i, vcov=vcovHC(i,type="HC0",cluster="countryname"))
summary(i2<-rdd_reg_lm(rdd_object= usa_f2,order=1)) 
coeftest(i2, vcov=vcovHC(i2,type="HC0",cluster="countryname"))

#nonparametric fuzzy
summary(rdbwselect_2014(y = manifesto_pop1$us_l, x = manifesto_pop1$dif_l1, c=0, bwselect="CCT"))
summary(rdrobust(y = manifesto_pop1$us_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, p=2, fuzzy = manifesto_pop1$treatment, h=1.383, b=3.036, covs=cbind(manifesto_pop1$state.d),cluster=manifesto_pop1$countryname, all=TRUE))
summary(rdrobust(y = manifesto_pop1$us_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, p=2, fuzzy = manifesto_pop1$treatment, h= 6.118, b=3.987, cluster=manifesto_pop1$countryname, all=TRUE))

# plot (Figure A1)
d1<-manifesto_pop1%>%dplyr::filter(dif_l1>-5, dif_l1<5)
pdf("Fig.A1b.pdf")
rdplot(x=d1$dif_l1, y=as.numeric(d1$us), covs=d1$state.d, c=0, p=2, x.label="", y.label="US Salience" , title="",y.lim = c(-0.1,0.3), x.lim = c(-5,5))
dev.off()
################################################################################
# Figure A4 plotting weighted estimates 
pdf("Fig.A4.pdf")
data <- data.frame(
  stringsAsFactors = FALSE,
  Models = c("CCT","CCT","CCT", 
             "MSE","MSE", "MSE",
             "CER", "CER","CER",
             "IK","IK","IK"),
  Estimate = c(0.465,  0.481, 0.481,
               0.430, 0.434, 0.434,
               0.480, 0.495, 0.495,
               0.430,0.457 ,0.457 ), 
  Method=c("Conventional", "Bias-Corrected", "Robust", 
           "Conventional", "Bias-Corrected", "Robust",
           "Conventional", "Bias-Corrected", "Robust",
           "Conventional", "Bias-Corrected", "Robust"),
  ord=c(1,2,3,4,5,6,7,8,9,10,11,12),
  CImax=c (0.754,0.770,0.797,
           0.568,0.572,0.620,
           0.762, 0.777, 0.792,
           0.534,0.561,0.612 ),
  CImin=c(0.176,0.192,0.166,
          0.292, 0.296, 0.248,
          0.198,0.213,0.198,
          0.325,0.352,0.301),
  se= c( 0.148,0.148,0.161,
         0.070, 0.070, 0.095,
         0.144, 0.144, 0.151,
         0.053 ,0.053,0.079))

data$ord <- factor(data$ord, levels = data$ord)

valXdodge = .25
ggplot(data = data, aes(x = Models, y = Estimate, group = Method)) +
  geom_point(position = position_dodge2(width = valXdodge), aes(shape=Method)) +
  geom_linerange(aes(ymin = CImin, ymax = CImax,linetype=Method),position = position_dodge(width = valXdodge))+labs(x="", y= "RD Estimates")+ ylim(0, 1)+
  coord_flip()+theme(legend.title = element_blank())

dev.off()



##########################################################################################

#-----------------------------------------------------------------------------
# DV= change of t and t-1 as Yi (first difference) 
#-----------------------------------------------------------------------------

#1. military position (Table A3)
manifesto_pop$military_change_l<-log(manifesto_pop$military_change+1)
manifesto_pop1<-manifesto_pop%>%dplyr::select(dif_l1, military_change_l,state.d, countryname,treatment)
manifesto_pop1<-manifesto_pop1[complete.cases(manifesto_pop1), ]
fakefuzzy <- rdd_data(x=manifesto_pop1$dif_l1,covar=cbind(manifesto_pop1$state.d), y=manifesto_pop1$military_change_l, z=manifesto_pop1$treatment,cutpoint=0)
fakefuzzy2 <- rdd_data(x=manifesto_pop1$dif_l1, y=manifesto_pop1$military_change_l,z=manifesto_pop1$treatment,cutpoint=0)

#parametric fuzzy
summary(f1<-rdd_reg_lm(rdd_object=fakefuzzy,covar.opt = list(strategy = c("include"), slope = c( "separate")), order=1))
coeftest(f1, vcov=vcovHC(f1,type="HC0",cluster="countryname"))
#w/o FEs 
summary(f2<-rdd_reg_lm(rdd_object=fakefuzzy2, order=1))
coeftest(f2, vcov=vcovHC(f2,type="HC0",cluster="countryname"))

#non-parametric fuzzy
summary(rdbwselect_2014(y = manifesto_pop1$military_change_l, x = manifesto_pop1$dif_l1, c=0, bwselect="CCT")) #IK
summary(rdrobust(y = manifesto_pop1$military_change_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, p=2, h=1.772,b=2.914, fuzzy = manifesto_pop1$treatment, covs=cbind(manifesto_pop1$state.d), cluster=manifesto_pop1$countryname, all=TRUE))#w/FE
summary(rdrobust(y = manifesto_pop1$military_change_l, x = manifesto_pop1$dif_l1, c=0, kernel = "tri",level = 95, p=2, h=1.772,b=2.914, fuzzy = manifesto_pop1$treatment,  cluster=manifesto_pop1$countryname, all=TRUE))#fuzzy w/o FE

#plot (Figure A2)
d1<-manifesto_pop1%>%dplyr::filter(dif_l1>-5, dif_l1<5)
pdf("Fig.A2.pdf")
rdplot(x=d1$dif_l1, y=as.numeric(d1$military_change_l), covs=d1$state.d, c=0, p=2, x.label="", y.label="Defence Position First Difference" , title="",y.lim = c(-1,1), x.lim = c(-5,5))
dev.off()

sink()



